Bond order solid of two-dimensional dipolar fermions 
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Recent experimental realization of dipolar Fermi gases near or below quantum degeneracy provides 
opportunity to engineer Hubbard-like models with long range interactions. Motivated by these 
experiments, we chart out the theoretical phase diagram of interacting dipolar fermions on the 
square lattice at zero temperature and half filling. We show that in addition to p-wave superfluid 
and charge density wave order, two new and exotic types of bond order emerge generically in dipolar 
fermion systems. These phases feature homogeneous density but periodic modulations of the kinetic 
hopping energy between nearest or next-nearest neighbors. Similar, but manifestly different, phases 
of two-dimensional correlated electrons have previously only been hypothesized and termed "density 
waves of nonzero angular momentum". Our results suggest that these phases can be constructed 
flexibly with dipolar fermions, using currently available experimental techniques. 

PACS numbers: 



Experimental demonstration of Bosc-Einstein conden- 
sation of atomic chromium [l[ and dysprosium 0], both 
of which have large magnetic dipole moments, ushers the 
ultra-cold dipolar gas to the arena of quantum emula- 
tion 0, A gas of the fermionic isotope of dyspro- 
sium 



acy 



161 Dy, has been cooled below quantum degener- 
. A high space-density gas of 40 K 87 Rb, fermionic 



molecules with electric dipole moments, has recently been 
produced near quantum degeneracy |6j and confined in 
optical lattice |7|. Such systems are expected to show 
a rich array of quantum phases arising from the long- 
range and anisotropic nature of dipole-dipole interaction 
[8l4lOl| . This uniquely distinguishes the dipolar Fermi gas 
from other Fermi systems, e.g. the 2D electron gas, the 
quantum fluid of 3 He, and Fermi gases of alkali atoms 
with short range interactions. Previous works on dipo- 
lar Fermi gases have investigated the anisotropic Fermi 
liquid properties 1(J U|, the pairing instability (l2- 16[. 
phases showing density modulation [UEij], as well as liq- 
uid crystal states III- 21 1. The possibility of supersolid 
phases [22| has also been discussed. 

For a 2D dipolar Fermi gas on a square lattice at half 
filling, with dipole moments perpendicular to the plane, 
one expects to find a checkerboard density modulation, 
known as the charge density wave (CDW, we follow the 
nomenclature even though atoms/molecules are charge 
neutral). When the dipole moments are aligned in the 
lattice plane the system becomes an anisotropic super- 
fluid and the attractive interaction binds fermions into 
Cooper pairs. The main question we address here is, how 
do different orders compete or cooperate as the dipole 
moments are turned from perpendicular to parallel ori- 
entation? 

We employ the functional renormalization group 
(FRG) technique ^MMk, 

along with self consistent mean 
field (SCMF) [26[ to obtain, for the first time, the zero- 



temperature phase diagram of dipolar fermions on a two 
dimensional lattice at half filling. The FRG takes an 
unbiased approach to treat all the instabilities of the 
Fermi surface, revealing the existence of two new and 
fascinating quantum phases: the p-wave bond order solid 
(BOSp); and the d-wave bond order solid (BOS^)- These 
bond order solids may be considered as 2D analogues of 
the "bond order wave" found in the ID extended Hub- 
bard model [i3-i3. 

We model single-component dipolar fermions on a two- 
dimensional square lattice with lattice constant ol by the 
Hamiltonian 
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where t represents the nearest neighbor hopping, ai is 
the fermion annihilation operator at the site i, rii = aja; 
is the number operator. The site index i represents a 
lattice site centered at r; = i x a^x + i v a\,y, where i x , 
i y are integers. The matrix elements of the dipole in- 
teraction in the two-particle Wannier basis are given by 
Vij = (ij\V dd \ij) = V d [l - 3(hj ■ d) 2 ]/{r l3 /a h f, where 
Tij = Ti — Tj and the dipoles are pointing in the same 
direction d. We assume an external electric or magnetic 
field F pointing in some general direction. Then the in- 
teraction energy of the dipole moment d with the field 
F is equal to — F • d, implying that the orientation of 
the dipole moments can be tuned by F. We label the 
direction of d by polar and azimuthal angles Op and (j>p 
respectively, as illustrated in the schematic of Fig. [If a). 

The interaction between dipoles can be attractive or 
repulsive depending on Op, cjj-p and r^. For example [re- 
fer to Fig. UK a)], if 4>f = 0, V y = Vdd(aLy) is always 
repulsive, while V x = Vdd(flLi) and V3 = Vdd{ah£ + o-lH) 
become negative for Of > i? c i ~ 35.26° and Op > t9 c2 = 
cos _1 (l/\/3) ~ 54.74° respectively. We shall show that 
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FIG. 1: (Color online) Dipolar fermions on square lattice, (a) Schematic of the dipolar fermions confined to a square opti- 
cal lattice potential. The induced dipole moment d points along the direction d = cosOfz + sin Of cos 4>fx + sin Of sin 0f?/. 
(b)Phase diagram obtained via FRG indicating four phases: p-wave bond order solid (BOS p ), d-wave bond order solid (BOSd), 
checkerboard charge density wave (cft-CDW), and p-wave BCS superfluid (BCS); left panel- phase diagram in the Op-Vd plane 
at <f>F = 0; right panel- phase diagram in the 0f-4>f plane at Vd = 0.5t. The phase boundary (solid line) is determined by the 
abrupt change in the symmetry of the eigenvector of the dominant instability (see Fig. [2}. The smooth crossover from cb-CDW 
and BOSd is indicated by a gradual change of the color shading, (c)-(e) Schematic of the bond or density modulation pattern 
for the BOSp, BOSd, and cfc-CDW phase respectively. 



these two critical points, i? cl and i9 c2 , roughly set the 
phase boundary between the checkerboard charge den- 
sity wave (c6-CDW), BOS p , and the Bardeen-Cooper- 
Schricffcr (BCS) superfluid phase, for the </>f = case. 

We now discuss the T = phase diagram at half 
filling. First, we analyze the weakly interacting limit, 
Vd < t, using FRG. In this approach, no assumptions 
about possible dominant orders are necessary. Rather, 
the method includes all processes near the Fermi surface 
of the non-interacting system via the generalized 4-point 
vertex function: Ui (ki, k 2 , k.3), where ki^ (ks^) are in- 
coming (outgoing) momenta and k4 = ki+k2 — kj. Here, 
£ is the rcnormalization group flow parameter that relates 
the energy cutoff A to the initial cutoff Ao (chosen to be 
At) via Ai = Age . Starting with the bare vertex Uo, 
progressively tracing out the high energy degrees of free- 
dom, a set of coupled intcgro-differcntial equations give 
the FRG flow for all the vertices. 

The renormalized vertex for specific channels of inter- 
est, e.g., 

C/ £ N E ST (ki,k 2 ) = C/,(k 1 ,k 2 ,k 1 + Q), " 
E/? cs (ki,k a ) = Ui(ki,-ki,k 2 ), 

are extracted by appropriately constraining the in- 
coming and out-going momenta. Here Q = (71", ±7r) is 
the nesting vector at half filling for the square lattice, 
and [/Nest is the game ag jjCDW of Ref Q xhe chan _ 

nel matrix with the largest divergent eigenvalue A cor- 
responds to the most dominant instability of the Fermi 
liquid. The corresponding eigenvector ip defined on the 
Fermi surface, indicates the symmetry of the incipient 
order parameter associated with the instability. 



(2) 



We perform the FRG analysis for a range of values of 
Vd, Of-, and cV producing a 3D phase diagram, visualized 
in Fig. [T{b) as slice cuts along two different planes. To 
capture and emphasize the key elements of the phase dia- 
gram, first we fix <j)-p = 0, generating a 2D phase diagram 
in the 9-p-Vd plane shown in the left panel of Fig. Hfb). 
Next we fix Vd = 0.5i instead, yielding the #f~</>f plane 
shown in the right panel of Fig. [ljb). 

The 6*p-Vd phase diagram shows the existence of three 
phases separated by two critical angles 9-p = 9\ and 62, 
with no appreciable dependence on Vd ■ For < 8p <6\, 
the nesting channel has the largest (most divergent) 
eigenvalue A. The corresponding eigenvector t/'nest, as 
illustrated in top panel of Fig. [2Ja), is almost constant 
with only small modulation along the Fermi surface. This 
implies the onset of CDW order with s-wave symmetry, 
identified as a checkerboard modulation of on-site den- 
sity, the cfr-CDW shown in Fig.QJe). The physical origin 
of this phase can be traced by observing that 9\ rts #i c , 
thus V x , V y , V3 > in this regime, allowing for a low 
energy configuration with density concentrated on the 
next-to-nearest neighbor sites, consistent with the per- 
fect nesting of the Fermi surface. For 62 < 9v < 90°, the 
BCS channel exhibiting a p-wave symmetry is the most 
diverging under FRG flow [see Fig. |2Ja)]. In real space, 
this corresponds to the onset of nearest neighbor pairing, 
(ai(Xi + x) = —{didi-a) generated by couplings V x and V3, 
both becoming attractive for #p > $2 ~ ^2c- The super- 
fluid phase here is the lattice analog of the p-wave BCS 
phase discussed previously for continuum dipolar Fermi 
gases @, [II Ell. 
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FIG. 2: (Color online) FRG results for Vd — 0.5t. The FRG is implemented numerically by discretizing the Fermi surface 
into 32 patches distributed at equally spaced angular points, (a) Top, middle and bottom panels represent FRG results for 
{0f,4>f) = (30°, 0), (42°, 0) and (70°, 0) respectively. Left column: the largest eigenvalue A of the NEST (dashed line) and BCS 
(solid line) channel. Right column: the corresponding eigenvector ip of the most diverging channel as function of £, the angle of 
the discrete k points on the Fermi surface defined by tan£ = k y /k x , plotted with square markers, (b) Top, middle and bottom 
panel represent FRG results for (#f,</>f) = (62°, 40°), (46°, 40°) and (38°, 40°) plotted using square markers. The fit is shown 
in solid line. As 9f is increased, ip smoothly changes from nodeless for 6f ;$ 46° to one with nodes for Of j5 46°. 



Finally the intermediate regime, 0\ < 6f < 02, is the 
most intriguing. The FRG predicts a leading instabil- 
ity in the nesting channel, similar to the cfo-CDW, but 
instead with a p-wave symmetry, '(/'nest (k) ~ x(k) = 
Xo sin k y , as shown in middle panel of Fig. HJa). This 
result suggests a broken symmetry phase, shown in 
Fig. [TJc), with periodic modulation of (ajdj+g — Xy) = 
-{a\a,i-y - Xy) = <5(-l) ia!+!y , where Xy is average of 
{a\a,i +y ) over all bonds. We observe that the nesting vec- 
tor Q is consistent with the checkerboard pattern of bond 
variable representing nearest-neighbor hopping. We re- 
fer to this broken symmetry phase as the p-wave bond 
order solid (BOS p ). Phases with similar, but manifestly 
different bond order patterns were conjectured by Nayak 
and referred to as p-density waves [3(J ■ 

The right panel of Fig. QJb), 0f~4>f phase diagram at 
fixed interaction strength, Vd = 0.5i, shows the three 
phases above for small values of </>f- However, as <pF 
is increased towards 45°, the BOS p region shrinks and 
eventually disappears beyond 4>f ~ 35°. Such change 
is due to the new features in the dipolar interactions 
for 0p close to 45°, where V x ~ V y , but the next-to- 
nearest neighbor interaction along x + y and x — y de- 
velop opposite sign. We find that for such large values 
of 4>f ~ 45°, the eigenvector can be fit very well by 
"0NEST (k) = a + /? [cos k x cos ky + sin k x sin k y ] , as seen 
in the right panel of Fig. [2]Jb). As Of is increased, the 
constant term a, which describes the density modula- 
tion of cfe-CDW order, is gradually reduced, while the 
magnitude of (3 increases. In the green shaded region 
in Fig. [Ub), a/f3 drops gradually from 1 to as 9p is 
increased toward the phase boundary to BCS. We re- 



fer to this region where the cos k x cos k y and sin k x sin k y 
components of t/'nest dominant as the d-wave bond or- 
der solid (BOSd). In this phase, the density and the 
nearest hopping {a\a i+x /y) are homogeneous. But the 
dipolar interaction induces an effective diagonal hopping, 
(a|ai_£ + y), a bond variable with amplitude proportional 
to /3 and spatial pattern shown schematically in Fig.QJd). 
BOSd found here differs from the cfcy-density wave con- 
jectured in Ref. [3p| ] - 

To firmly pin down the nature of the phases, we com- 
plement the FRG analysis with SCMF theory (see Ref. 
[26( ) on a square lattice of finite size L x L with period 
boundary condition by defining the normal and pair den- 
sity matrices pij = (ajoti) and n%ij = (a,ia,j) respectively. 
The corresponding mean fields are then given by Xji = 
-J2ki(j k \ v d<i\ li )Pik and Aij = -| X) fc j(ij|Vdd|fcO m tt- 
The dipole interaction is retained up to a distance of 
IOol. We search for the ground state iteratively by start- 
ing with an initial guess for p and m, until desired con- 
vergence is reached. The phase boundaries are obtained 
by comparing the thermodynamic potential for various 
converged solutions (see Supplementary Material). The 
chemical potential is tuned to maintain half filling. And 
the lattice size L > 20a^ is varied to check the results do 
not depend on the choice of L. 

The SCMF phase diagram for 4>p = 0, shown in Fig. [31 
confirms the existence and interpretation of the three 
phases found in the FRG analysis. The phase boundaries 
are in qualitative agreement with those from FRG. SCMF 
for non-zero 4>f also identifies the BOS^ as a phase with 
the bond modulation pattern illustrated in Fig.[IJd). We 
caution that the SCMF phase diagram is only suggestive. 
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FIG. 3: (Color online) SCMF phase diagram. Shown on the 
left are representatives of the on-site density pa, the nearest 
neighbor hopping p^ (with j = i + x or j — i + y), or the 
pairing gap Ajj corresponding to the four phases at Va = 0.51 
Lattice size is 32 x 32. 



For example, SCMF predicts an additional striped den- 
sity wave phase, the st-CDW, which is not expected to 
survive at Vd <C t. This illustrates that SCMF is insuf- 
ficient to describe competing orders as opposed to FRG. 
The possibility of st-CDW and collapse instability be- 
yond the weak coupling regime is further discussed in 
the supplementary material. 

We now provide some intuitive understanding of the 
bond order phases by considering a simplified mean field 
version of Eq. (UJ , keeping only the nearest neighbor in- 
teractions V x and Vy. The mean field decoupling of the 
interaction term gives — n,rij ~ a\aja^ai — > />, ,a\a, + 
h.c. — \pij\ 2 . The modulation of the bond variable, pij = 



(ajoj), in the BOS p phase at 4>f = has the form show 



in Fig. [He), p l . i ± x = Xx, Pi.i± y = Xy±S- The mean field 
Hamiltonian can be written as Hr = — 2^ k Xkfri c ak + 
h.c, up to a constant term. Here and 6k are fermion 
annihilation operators defined separately on two sub- 
lattices related by the lattice translation vector a^x, and 
Xk = (t + VxXx) cos k x + (t + VyXy) cos k y - iV y S sin k y . 
The ground state energy per unit cell is then given by 
£gs = ~2(Xx +Xy)(t + V x + V y ) - 2V y S 2 , clearly indicat- 
ing that finite bond modulation S is energetically favored 
for positive V y . The 4>f = 90° situation is identical, only 
with x and y axis interchanged, and hence a 90° rotated 
bond pattern. Thus, the BOS<j phase, with checkerboard 
pattern of next-to- nearest bonds near <j)p = 45°, natu- 
rally connects the two BOS p phases on either side. 

The bond modulation 5, the energy gap, and the tran- 
sition temperature T c of the BOS p phase increase with Vd 
for weak coupling. Exact diagonalization of Eq. (JXJ) on a 
2x8 and 4x4 cluster with periodic boundary conditions 
shows that the optimal place to observe the BOS„ is at 



intermediate interaction and tilt angle, e.g. Vd ~ 2.5t 
and (Op, 4>f) = (45°, 0°), where the energy gap, and thus 
T c , is maximal. Mean field theory estimates an optimal 
T c - 0.23t, or about 0.05E F for half filling, which is not 
too far from the temperature achieved in Dy experiment, 
T - 0.25£ F [!■ The BOS d on the other hand is most sta- 
ble in the vicinity of 4>f = 45° for 9f ~ 60°. The charac- 
teristic density modulation of the cfr-CDW and st-CDW 
phase uniquely distinguishes them from the other phases 
and may be detected via in-situ density imaging. The 
BCS phase can be detected via pair correlation measure- 
ments using noise spectroscopy [3l|. Finally the BOS^ 
phase may be distinguished from BOS p by probing the 
d-wave symmetry via the pump-probe scheme discussed 
in Rcf. [32j. Finally, in the presence of a trap potential, 
the insulating plateau at half filling will be surrounded 
by metallic regions. The approaches outlined here can 
be employed to study dipolar Fermi gas away from half- 
filling. 
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FUNCTIONAL RENORMALIZATION GROUP 

Functional renormalization group (FRG) is a powerful 
approach to correlated fcrmion systems. In recent years, 
it has been applied successfully to a wide of range of 
problems to yield decisive insights impossible to obtain 
from other approaches [l| . Our implementation of FRG 
follows Zanchi and Schulz 0], which can be viewed as a 
generalization of earlier renormalization group approach 
to interacting fermions, e.g., by Shankar |3j and Polchin- 
ski Q ■ Our implementation has been employed to study, 
for example, the quantum phases of Bose- Fermi mixtures 

1- 

Here we provide further technical details of our nu- 
merical FRG calculation. The Fermi surface of the 
non-interacting system is a nested square for half fill- 
ing. It is discretized into TV = 32 patches with mo- 
menta ki, i G {1,2,.. TV}. The 4-point interaction ver- 
tex as a function of three momenta U( (ki, k2, ka) is 
represented as a N 3 x N 3 matrix. We start with the 
bare vertices defined on the discretized Fermi surface, 
L/e =0 (ki,k.,-,ki) = Vdd(ki-ki), where V d d(q) is the dipole 
interaction in momentum space, i.e., the lattice Fourier 
transform of Vdd(i"ij) given in the main text, with proper 
anti-symmetrization required by Fermi statistics. The 
momentum dependence of the dipolar interaction is fully 
taken into account. At each FRG step £, the renormalized 
vertices Ui(ki, k,-, k;) are calculated at one- loop level, 
truncating the effective action at the six fermion terms 
and neglecting self-energy correction The renormal- 
ized vertices for particular channels, e.g. Uf EST and 
Uf cs , are then extracted by appropriately choosing the 
in-coming and out-going momenta. More technical de- 
tails on the flow equations can be found in Ref. Q. We 
emphasize that the traditional single-channel (in cither 
particle-particle or particle- hole channel) RG, or Random 
Phase Approximation, is insufficient to reveal complex 
order such as BOS p and BOS^ as discussed here. 

Within our implementation of FRG, the shape of the 
Fermi surface is fixed. Self-energy corrections in princi- 
ple can distort the Fermi surface, potentially making the 
nesting at Q = (ir, it) less perfect or even favoring other 
nesting such as (0, it). These are higher order effects that 
can be safely neglected in the weak coupling limit. Sim- 
ilar arguments were made in many FRG studies e.g. on 
the Hubbard model Q • The FRG phase diagram is most 



reliable and trustworthy in the weak coupling regime, 
Vd « t. The phase diagram is expected to be modified 
at stronger interactions. 

SELF-CONSISTENT MEAN FIELD THEORY 

The self-consistent mean field theory adopted here fol- 
lows Ripka et al Q. It allows coexistence of multiple 
orders (in particular bond and current order) and long 
range interactions. We define the real space generalized 
density matrix on a square lattice of size L x L with 
periodic boundary conditions, 



P 

m* 



rn 
1 p* 



(1) 



The matrices p and m are of dimension L 2 x I? and 
defined as pij = (ajaj) and rriij = (<XiOj), and represent 
normal and anomalous averages, respectively. 1 is the 
L 2 x L 2 identity matrix. 1Z is thus of dimension (2L 2 ) x 
(2L 2 ). Note that include both on-site density pa and 
nearest neighbor hopping, e.g., pi.i+x- The mean field 
Hamiltonian of Eq. (1) in the main text is 



H m f = 



-A* 



A 

lUki^kimik and Xji 



(2) 



with the mean fields Aji 
— ^luVjkUPik- t is the tunneling matrix, i.e., Uj has 
the value t if the sites i and j are connected by a bond. 
We search for the ground state iteratively by starting 
with an initial guess for p and m. Each iteration con- 
sists of three steps: (i) construct H m f using the defini- 
tions for x an( i A above, (ii) diagonalize H m f to ob- 
tain the eigenvectors W v with negative energy eigen- 
values E„, U mf W v = -\E U \W V , and (iii) update 
and m by constructing 1Z via 7Z = J2„ W V W V ' 
The iteration is repeated until desired convergence is 
reached. The phase boundary is determined by com- 
paring the thermodynamic potential of different con- 
verged solutions, T = Tr [(— t — pi + x)P ~ Am*] — 
3 Hijki V~ijki(2pkiPij - m^mik). The chemical potential 
is tuned to maintain half filling. And the lattice size L is 
varied to check the results do not depend on the choice 
of L. 

We caution that the SCMF phase diagram is only sug- 
gestive. It shows candidate phases that are energetically 
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competitive on the mean field level. On one hand, it 
can take into account the self-energy corrections (on the 
Hartree-Fock-Bogoliubov level) which become increas- 
ingly important for larger values of Vd/t. It is assuring 
to see that BOS p persist. On the other hand, it is not 
always guaranteed to provide accurate predictions espe- 
cially for intermediate dipole tilt angles where different 
orders compete. For example, it misleadingly predicts st- 
CDW in the weak coupling limit Va <C t. In the striped 
phase, the on-site density is modulated with period 2a^ 
in the y direction while remaining translationally invari- 
ant along x. We find its mean-field energy to be very close 
to the BOSp phase. However, it is important to note st- 
CDW is ruled out in the weak coupling limit by FRG and 
on physical grounds that the Fermi surface deformation 
required to favor the (0, 7r) nesting vector becomes neg- 
ligible at weak coupling. In other words, we expect that 
fluctuations beyond mean field will shrink the region of 
st-CDW so it cannot survive to Vd « t. 

Finally, we comment on the possibilities of st-CDW 
and collapse instability at finite Vd/t, beyond weak cou- 



pling. Neither FRG or SCMF can unambiguously settle 
these issues. And we have to wait for future large scale 
exact diagonalization or quantum Monte Carlo studies. 
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